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Abstract 

We investigate a method to compute a finite set of preliminary orbits 
for solar system bodies using the first integrals of the Kepler problem. This 
method is thought for the applications to the modern sets of astrometric 
observations, where often the information contained in the observations 
allows only to compute, by interpolation, two angular positions of the 
observed body and their time derivatives at a given epoch; we call this 
set of data attributable. Given two attributables of the same body at two 
different epochs we can use the energy and angular momentum integrals 
of the two-body problem to write a system of polynomial equations for 
the topocentric distance and the radial velocity at the two epochs. We 
define two different algorithms for the computation of the solutions, based 
on different ways to perform elimination of variables and obtain a univari- 
ate polynomial. Moreover we use the redundancy of the data to test the 
hypothesis that two attributables belong to the same body (linkage prob- 
lem). It is also possible to compute a covariance matrix, describing the 
uncertainty of the preliminary orbits which results from the observation 
error statistics. The performance of this method has been investigated by 
using a large set of simulated observations of the Pan-STARRS project. 



1 Introduction 

With the new observational techniques of the next generation surveys, like Pan- 
STARRS and LSST0 the number of moving objects detected in each night of 
observations is expected to increase by two orders of magnitude with respect to 
the current surveys. To deal with this huge amount of data the interest in the 
study of orbit determination methods has been renewed, both from the theoret- 
ical and the computational point of view. The classical methods of preliminary 

*Giovanni F. Gronchi: Dipartimento di Matcmatica, Universita di Pisa, Largo B. Pon- 
tecorvo, 5, Pisa, Italy gronchiOdm.unipi.it 

^Linda Dimare: Dipartimento di Matematica, Universita di Roma 'La Sapienza', P.le Aldo 
Moro, 2, Roma, Italy dimare@mat.uniromal.it 

* Andrea Milani: Dipartimento di Matematica, Universita di Pisa, Largo B. Pontecorvo, 5, 
Pisa, Italy milani@dm.unipi.it 

1 see the web pages http://pan-starrs.ifa.hawaii.edu, and http://www.lsst.org 



orbit determination by Laplace [5] and Gauss [3], that have been often revis- 
ited in the last two centuries [18], [9J, [TO], are based on the knowledge of at 
least three observations of a solar system body in three different nights. Both 
Laplace's and Gauss' method may produce more than one preliminary orbit for 
the same object: a detailed analysis of the occurrence of multiple solutions is 
in [B]. The determination of a preliminary orbit is followed by the differential 
corrections [TJ , an iterative method to obtain the minimum of a target function, 
that improves the orbit in the sense of the least squares fit of the residuals: this 
sequence of operations was already proposed in [3]. 

The data of the current surveys generally do not provide a single observation 
for an object in an observing night: in fact the moving objects are distinguished 
from fixed stars by detecting them a few times in the same night: the sequence 
of observations usually gives a too short arc on the celestial sphere [14], such 
that the data are not enough to compute an orbit for that body. As the number 
of detected objects per night is very large, it is difficult to decide whether two 
sequences of observations made in different nights belong to the same object: 
this gives rise to the problem of linkage of two short arcs of observations. 

The information contained in a short arc of observations can be used to define 
an attributable [TTJ, consisting of the angular position and velocity of the body 
on the celestial sphere at a given time; the topocentric distance and the radial 
velocity at that time are unknown. Therefore two short arcs of observations 
belonging to the same object provide us with 8 scalar data, from which we can 
try to compute an orbit. 

In 1977 Taff and Hall proposed to use the angular momentum and the en- 
ergy integrals to perform orbit determination starting from a data set that cor- 
responds to two attributables of the same observed body (see [TO] . [20]). They 
noticed that the problem can be written in an algebraic form but, since the 
total degree is high, they suggested to use a Newton- Raphson method to solve 
the problem. This approach deals with the solutions only locally, and there are 
alternative possible solutions that can be lost. In this paper we shall start from 
the same first integrals of the Kepler problem, but we shall exploit the algebraic 
character of the problem, keeping in this way a global control on the solutions. 
In particular we shall present two different methods to solve by elimination 
the polynomial system corresponding to this problem, and to compute all the 
related preliminary orbits defined by the two attributables. 

Since an orbit is defined by 6 scalar data, the available information is re- 
dundant, and we can use this to set compatibility conditions for the solutions 
(see |H])), that should be fulfilled if the two attributables belong to the same 
solar system object. The unavoidable errors in the observations affect also the 
computation of the attributables. Given a covariance matrix for the two at- 
tributables, expressing their uncertainty, we can use this to compute the value 
of an identification norm, based on the compatibility conditions, to decide if the 
attributables may be related to the same body (i.e. if the linkage is successful) 
and to choose among possible alternative solutions. 

The plan of the paper is the following: after introducing some notation 
related to attributables in Section [21 we explain in Section [3] the orbit determi- 
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nation method and derive the bivariate polynomial system whose roots give us 
the topocentric distances of the observed object at the two epochs. In Section|4] 
we explain the two algorithms that we propose to search for these roots and the 
effective computation of the orbits. In Section[S]we deal with the uncertainty of 
the data, introduce the identification norm and provide the covariance matrices 
of the preliminary orbits. Some numerical experiments are presented in the last 
two sections: in Section [6] we show the results of a test case, illustrating the 
computation of the preliminary orbits for a numbered asteroid whose orbit is 
well known, while in Section [7] we investigate the performance of the method 
for a large database of simulated observations. 

2 Attributables 

Let (p,a,5) 6 R + x [— tt, tt) x (— 7r/2,7r/2) be spherical coordinates for the 
topocentric position of a solar system body. The angular coordinates (a, 5) 
are defined by a topocentric coordinate system that can be arbitrarily selected. 
Usually, in the applications, a is the right ascension and S the declination with 
respect to an equatorial coordinate system (e.g., J2000). 

Given a short arc of observations of a celestial body (i,, aj. Si), for z = 1 ... to 
with to > 2, it is often possible to compute an attributable^ that is a vector 

A = (a, 5, a, 5) 6 [-tt, tt) x (-tt/2, tt/2) x R 2 , 

representing the angular position and velocity of the body at a mean time t in 
the selected coordinates (see [H]). Usually we choose t as the mean ti)/m. 
The attributable is computed by a polynomial fit, typically linear or quadratic, 
and the observations used in the computation need to be made by the same 
observatory. If the observations are enough, i.e. m > 2 for a linear fit, to > 3 for 
a quadratic one, then we can compute also a covariance matrix T^, representing 
the uncertainty of the attributable. Note that the topocentric distances pi at 
times tj are completely unknown. 

We introduce the heliocentric position and velocity of the body at time i 

r = q + pp, r = q + pp + p{p a a + psS) , (1) 

with p, p the topocentric distance and the radial velocity, and with p, p a , p$ 
the observation direction and its partial derivatives with respect to a and S. 

The vectors q, q represent the heliocentric position and velocity of the ob- 
server on the Earth. The observer position is know as a function of time, but 
for consistency, if the attributable is computed by a fit to polynomials with 
low degree, the values q(t),q(t) need to be computed by the same interpola- 
tion. Therefore we make a quadratic fit with the actual geocentric positions 
q(tj) — q©(£j) at the times of the individual observations (q® is the heliocentric 
position of the Earth center) to obtain the interpolating function q b s (£); then 

2 The name refers to the possibility of attributing the observations of the short arc to an 
already known orbit. 
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we take q(i) = q®(t) + q 6s(*) and q(t) = q©(i) + q 6s(*)- This method was 
suggested by Poincare in |18) . and it is important to obtain preliminary orbits 
of better quality, see [T5] . 

In rectangular coordinates we have 

p = (cos a cos S, sin a cos 5, sin 5) , 
p a = {— sin a cos 5, cos a cos S, 0) , 
ps = (— cos a sin <5, — sin a sin <5, cos 5) . 

These vectors form an orthogonal system, in particular 

|p| = |pi| = l, |/5 Q |=COS(5, p ■ p a = p ■ p S = Pa ' f>6 = , 

where the dot indicates the Euclidean scalar product and | • | the corresponding 
norm. 

For later reference we introduce the orthonormal basis {p, v,n} adapted to 
the apparent path p — p(t) of the observed body on the celestial sphere: we 
define the unit vector v by the relation 

d , 

where 77 = \J a 2 cos 2 8 + S 2 and is called proper motion. Moreover we set n = 
p x v. 



3 Linkage by the two-body integrals 

Given two attributables A\, A2 at different epochs ti, ?2, in the hypothesis that 
they belong to the same observed body, we write down polynomial equations 
for the topocentric distance and radial velocity of the body at the two epochs 
by using the angular momentum and the energy integrals. 

3.1 Angular momentum and Energy 

For a given attributable A the angular momentum vector (per unit mass) can 
be written as a polynomial function of the radial distance and velocity p, p: 

c(p, p) = r x r = Dp + Ep 2 + Fp + G , 

where 

D = q x p, 

E = ap x pa + 5p x ps = i]h , 
F = aq x p a + <5q x p s + p x q , 
G = q x q, 

depend only on the attributable A and on the motion of the observer q, q at 
the time t of the attributable. For the given A we can also write the two-body 
energy as a function of p, p, as in |12j 

2fc 2 

2£ (p, p) = p 2 + cip + c 2 p 2 + c 3 p + c 4 , = , 

VP + C 5P + CO 
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where k is Gauss' constant and 

co = |q| 2 , ci = 2qp, c 2 = 77 2 , 

c 3 = 2(d q • p Q + <5 q- p 5 ) , c 4 = |q| 2 , c 5 = 2qp, 

depend only on A, q, q. 

3.2 Equating the integrals 

Now we take two attributables A\ = (ai, &i, 5i), A2 = («2, ^2, «2, ^2) at 
epochs ti, i 2 ; we shall use the notation of Section ETT1 with index 1 or 2 referring 
to the epoch. If Ai, A2 correspond to the same physical object, then the angular 
momentum vectors at the two epochs must coincide: 

D 1 p 1 -D 2 p 2 =J(p 1 ,p 2 ), (2) 

where 

J(pi, Pa) = E2P2 - Eip 2 + F 2( o 2 - Fipi + G 2 - d . 

Relation @ is a system of three equations in the four unknowns pi,pi,p2,p2i 
with constraints 

pi > , p 2 > . 

By scalar multiplication of @ with Di x D 2 we eliminate the variables p\,p2 
and obtain the equation 

Di xD 2 -J(pi,pa) = 0. (3) 

The left hand side in ([3]) is a quadratic form in the variables p±, P2] we write it 

as 

q(Pl,P2) = <720Pl + 910/01 + 902P2 + 901/02 + <7oo , (4) 



with 



q 20 = -Ei Di x D 2 , 902 = E 2 • Di x D 2 , 

q w = -Fx Di x D 2 , q 01 = F 2 • Di x D 2 , 



9oo = (G 2 - Gi) • Di xD 2 . 

Equation (|4]) defines a conic section in the (pi,p 2 ) plane, with symmetry axes 
parallel to the coordinate axes. Since the directions of Ei,E 2 correspond to 
fii, n 2 , for |f 2 — ti| small enough the angle between these two directions is small 
and the coefficients 9 2 o,9o2 have opposite signs, thus in this case (U]) defines a 
hyperbola. 

We can compute the radial velocities p\ , p 2 by vector multiplication of ^ 
with Di and D 2 , projecting on the direction of Di x D 2 : 



. (J x D 2 ) ■ (Di x D 2 ) (J x Di) • (Di x D 2 ) 
P1(P1 ' P2) = | Dl xD 2 P ' = | Dl xD 2 |2 



(5) 
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For the given A\ , A2 we can also equate the corresponding two-body energies 
£1,82- We use the expressions of px{px, P2), f>2{px, P2) above and substitute them 
into Z\ = £2, thus we obtain 

•T^Oi,^)- l k = =-F 2 (pi,p 2 ) - l k , (6) 

for some polynomial functions J-i{pi,P2), Fzipi, P2), Gx(px), ^2(^2) with degrees 
deg(J r i) = deg(J r 2 ) = 4 and deg(^i) = deg(£ 2 ) = 2. By squaring we have 



(J 7 ! - ^ 2 ) 2 eie 2 - 4fc 4 (^ + &) = -8k 4 ^g 2 ■ (7) 
Squaring again we obtain the polynomial equation 

p(pup 2 ) = [(J 7 ! - TifQxGi - 4k\G 1 +g 2 )] 2 - 64fc 8 g 1 £ 2 = 0, (8) 

with total degree 24. Some spurious solutions may have been added as a result 
of squaring expressions with unknown sign. 

Note that, if the observations were made from the center of the Earth, 
Gi {i — 1,2) would be the angular momentum of the Earth at epochs t\ 1 i 2 , 
thus Gi = G 2 and goo — 0. With this simplifying assumption p\ = p 2 = 
is a solution of the system g(pi,p 2 ) = p(pi,p 2 ) = that corresponds to the 
Earth center and therefore is not acceptable. This solution also appears in the 
geocentric version of the method of Laplace for a preliminary orbit from 3 obser- 
vations. Actually we use topocentric observations, for which the zero solution 
is replaced by one with both pi and p 2 very small. 



3.3 Degenerate cases 

The quadratic form (j4]) degenerates into a linear function when 

Ei • Di x D 2 = E 2 • Di x D 2 = . 

A simple computation shows that 

Ei • Di x D 2 = ?/i(ni • qi)(pi x p 2 ■ q 2 ) , 
E 2 • Di x D 2 = ?7 2 (n 2 • q 2 )(pi x p 2 ■ qi) , 

thus, assuming that the proper motions i]2 do not vanish and setting ni2 = 
Pi x p2, the degeneration occurs when either n 12 vanishes (CO) or at least one 
the following relations holds: 

fii • qi = n 2 ■ q 2 = , (CI) 

n i2 • qi = n i2 • q 2 = 0, (C2) 

fii • qi = ni2 • qi = , (C3) 

n 2 • q 2 — ni2 • q 2 = . (C4) 

The interpretation of these conditions is the following: (CO) means that pi,p 2 
point to either exactly the same or exactly the opposite direction in the sky; (CI) 
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means that both sets {pi,vi,qi} and {p 2 ,v 2 ,q 2 } are constituted of coplanar 
vectors; (C2) says that pi,p 2 ,qi,q 2 are coplanar. Let us discuss condition 
(C3): fii • qi = means that qi,pi,vi are coplanar and 1112 • qi = means 
that pi,p 2 ,qi are coplanar as well. If Di ^ we obtain that the four vectors 
qi, px, vi, p2 all lie in the same plane. In particular (C3) implies that p 2 belongs 
to the great circle defined by the intersection of the plane generated by px,V\ 
with the celestial sphere. This degeneration condition can be compared with 
the failure condition of the classical orbit determination methods with three 
observations by Gauss and Laplace [T7], due to vanishing of the curvature in 
the apparent path of the observed body on the celestial sphere. The discussion of 
condition (C4) is similar to the previous one and corresponds to the coplanarity 
of q 2 ,p 2 ,v 2 ,pi. 

4 Computation of the solutions 

In this section we introduce two different methods to search for the solutions of 
the semi-algebraic problem 

1 q(puP2)=0 ' Pl ' P2> ° (9) 

for the polynomials p, q introduced in ©, ((H) respectively. Moreover we explain 
the full procedure for the computation of the preliminary orbits and introduce 
compatibility conditions to decide whether the attributables used to define the 
problem are related to the same solar system body. 

4.1 Computation of the resultant via DFT 

The first method consists in writing the resultant (see [3]) of p and q with 
respect to one variable, say p\. In this way we find a univariate polynomial in 
the p2 variable whose real positive roots are the only possible /9 2 -components of 
a solution of ([9]) . By grouping the monomials with the same power of p\ we can 
write 

20 

p{pi,P2) = ^2 a,j(p2) Pi , where (10) 

3=0 

(20 for j = ... 4 

deg(Oj-) = < 24-(j + l) forj = 2fc-l with k > 3 
[ 24 - j for j = 2k with k > 3 

and 

q(pi,P2) = h p\ + h pi + b (p 2 ) (11) 

for some univariate polynomial coefficients a*, bj depending on p 2 (actually 61, 6 2 
are constant). We consider the resultant i?es(p 2 ) of p, q with respect to p\\ it is 
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generically a degree 48 polynomial denned as the determinant of the Sylvester 
matrix 



«20 
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(12) 



a ai : : : b &i 
\0a 00 6 / 

The positive real roots of Res(p2) are the only possible values of p 2 for a solution 
(pi,/92) of ([§]). We could use the resultant method to eliminate the variable P2 
by a different grouping of the terms of p, q: 

20 

P{PI,P2) =^4(pi) p 3 2 , q{pi,Pi) = b' 2 p\ + b[ p 2 + b Q (pi) , 



where the degrees of a'j are described by the same rules as for aj . 

Apart from non-real and non-positive solutions, we shall see that there are 
additional different reasons to discard some pairs of solutions of ©, thus we 
expect that the number of acceptable ones is not large. We use a scheme similar 
to [5] to compute the coefficients of the resultant Res(p2)- 



1) evaluate a,i(p2), bj{p2) at the 64-th roots of unit w& 



= 27ri 7 



,k = 0...63, 



by a DFT (Discrete Fourier Transform) algorithm; 

2) compute the determinant of the 64 Sylvester matrices; by relation 

det(S(p 2 )| P2 =cuJ = (detS(p 2 )) \ P2 =u k 
we have the values of Res(p 2 ) at the 64-th roots of unit; 

3) apply an IDFT (Inverse Discrete Fourier Transform) algorithm to obtain 
the coefficients of Res(p2) from its evaluations. 

The use of the DFT and IDFT allows us to interpolate the resultant Res(p 2 ) 
in an efficient way. The use of numerical evaluations, e.g. at the roots of unit, 
avoids the difficulty of writing a very long symbolic expression for the resultant, 
that could be cumbersome to be managed by a programming language compiler. 

The complete set of complex roots of Res(p2), with an error bound for each 
of them, are computed using the algorithm described in [2], which is based on 
simultaneous iterations. Let P2{k), k — 1 . . .n < 48, be the subset of the real 
and positive roots of Res(p2). Then for each k we perform the sequence of 
operations below: 

4) solve the equation q(pi, p 2 (^)) = and compute the two possible values 
Pi(k, 1), pi(k, 2) for pi, discarding negative solutions. Then define pi(k) 
equal to either pi(k, 1) or pi(k, 2), selecting the one that gives the smaller 
value of \p(pi,p 2 (k))\; 



8 



5) discard spurious solutions, resulting from the squaring used to reduce the 
energy equality to the polynomial equation (|8|). The spurious solutions 
are the solutions of ^ that do not satisfy either (J7]) or ©; 

6) compute the corresponding values of f>i(k), f>2(k) by ([5]) and obtain a pair 
of orbits defined by the sets (ct^, <Jj, a,, pi, pi) of attributable elements^ 
for i = 1,2; 

7) change from attributable elements to Cartesian heliocentric coordinates by 
relation = pi(k)pi + for i = 1,2, and by the corresponding formula 
for fj. Note that the observer position is not the actual q(ti), but is 
obtained by interpolation as proposed by Poincare (see Section [2]). Then a 
standard coordinate change allows us to obtain the related pairs of orbital 
elements: we shall use Keplerian elements (a, e, J, Q,u),£), where I is the 
mean anomaljQ- The epochs of the orbits are i\(k),ii{k), corrected by 
aberration due to the finite velocity of the light c: ti(k) — U ~ pi(k)/c for 
i = 1,2. 

We have implemented this algorithm in FORTRAN 90 using quadruple pre- 
cision for part of these computations, in particular the ones related to DFT and 
IDFT. This feature appeared necessary to obtain reliable results starting from 
our first numerical experiments. 



4.2 Normal form of the problem 

Another method to compute the solutions of ([§]) is based on a coordinate change 
to variables (£i, £2)1 that allows to perform easily the elimination of either £1 or 
£2- Let us set 

20 

p(p%,p2) = ^2 Vi,ip\pi ■ 

i,j=0 

First we consider the affine transformation to intermediate variables (Ci, C2) 



T ' ' P2 J ~" V C2 J \ o- 2 1 p 2 -T 2 
where, to eliminate the linear terms in Q, we set 

9l,0 dcf 90,1 dcf a 

0-1T1 = = a , ct 2 t 2 = = p , 

2^2,0 2q ,2 



3 The attributable elements are the same as spherical polar coordinates with their time 
derivatives: the coordinates are just reordered in such a way that the first four elements form 
the attributable, hence the name. 

4 Any other set of orbital elements in which the first four are defined by the two-body 
energy and angular momentum can be used, e.g. cometary elements (p^, e, /, Q, ui, t p ) where 
Pd is the perihelion distance and t p the time of perihelion passage: this set would allow to 
handle also parabolic and hyperbolic orbits. 
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so that 



r 2 , 90,2O 2 2 K 
^l" 1 "" l2 t '2 + 



92,0^! 92,00"! J 

If we set, for an arbitrary a 2 E M, 

o-i = 702 , 7 = 

we obtain 

9 o T-^Ci, C2) = 92,oo? [Ci 2 - C 2 2 - 2c,] , 



with k = 9o,o — 



qIo «o,i 



4<?2,0 4g ,2 



go, 2 

92,0 



with c* 



2 92,00-1 



2 • 



We already observed that, for \t 2 — h\ small enough, q 02 and 920 have opposite 
signs, hence in this case the variable change T is real. However, in general, we 
have to consider T as a transformation of the complex domain C 2 . We also have 



P°r- 1 (Ci,c 2 )= ^2 Mid, 

i,j=0 



where 



20 20 



a h - i l3 k - j 



h—i k—j 

and a, /3, 7 depend only on the coefficients of q(pi,p 2 ). 

Now we apply a rotation of angle 7r/4 to pass to the (£1,^2) variables: 

n : 



Ci 

C2 



£1 

6 



cos(f) -sin(f) 
sin(f) cos(f) 



Ci 

C2 



We have 



9 oX" 1 oft" 1 (£1,6) 
poT-'oTl- 1 ^,^) 



2 92, oo 2 [£16 - c*] 

24 

E pSj&* 



where 
and 

Pm,n 







if i + j < 24 
if i + j > 24 



(13) 



E ?M E ( \ 

h+k=m i+j—n ^ 

min{m,20} n 

Ph,m-h 

/i=max{m-20,0} i=0 



j / 2(' l + fc )/ 2 



E^-^-^E 



/i=0 



i=0 



m — h 
n — i 



(-iy 



2'm/2 



(14) 



m-A\ (-1)"- 
n — i ) 2 m / 2 



10 



The last equality is obtained taking into account that pn.m-h = for h > 
20 or to — h > 20. Using the relation £1^2 = c* we can consider in place of 
poT~ x o TZ~ 1 (!;i,t;2) the polynomial 



24 24 24 



~* ( t t \ \ ^ * k t-h — k 1 \ ^ * ft 1 \ ^ * h t-k — h 

p (si, 6) = 2^ ^,fe c *^i + z^m> c * + 2^ PhM c ^2 = 

h, k = fr — h, fc = 

ft. > fc h < k 



24 



El E pU<£- 3 k?+£( E pU c *"'KI+EpU 



1 \ fe, fe = / \ fe,fe = o / 

h — k = j k — h = j 



with 



24 



A 3 = E PlA~ j = Y,pU-A'^ i = l..-24, 



/i, fc = /l— J 

t 24 



4>(&) - s 24 e 2 24 + . . . + B± 6 + B , 

24 24 

^= E = Y,p*k-iA~ j > i = o...24 



h — h = 3 



We consider the algebraic problem in normal form 

p*(£i,6) = o 

ei6 - = • 



(15) 



In this case we have to consider all the solutions of (fT5j) , not only the ones with 
real and positive components. 

If c* = 0, then the solutions (£1, £2) of (T5"j) are of the form (£i(fc), 0) or (0, &{k)), 
where £i(fc),6(fc), k = 1 • • ■ 24 are the roots of A 24 ,^f i + ... + Aifi + S and 
-B24 £| 4 + • ■ • + Bi£ 2 + B respectively. 

If c* 7^ 0, using the relation £1^2 = c* we can eliminate one variable, say £1, 
from p* . Thus we obtain the univariate polynomial 

ls ' A 24 _ fc cf- fe , < fc < 23 



P(6) = E^, with p fc = \ 

k=0 



^fe-24 , 24 < fc < 48 



We compute all the complex roots £2(^)1 k = 1 ... 48 of p(£ 2 ) by the algorithm 
in [2] ; then for each k we define the other component of the solution by 

£i(*0 



6(fc) 



Given all the complex solutions of (|15[) we compute the corresponding points in 
the (pi,p2) plane by 

(pi(fc),p 2 (fc)) = T- 1 on- 1 (a(fc),6(fc)) , fc = i...48, 
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discarding the ones with non-real or non-positive components. At this point the 
preliminary orbits can be computed following the same steps 5), 6), 7) of the 
algorithm explained in Subsection 14.11 

From a few experiments performed this method seems to require more than 
quadruple precision because of the complicated formulae defining the transfor- 
mation used to obtain the normal form (|15p . Thus the advantage in the simple 
elimination of the variable £i must be balanced with the introduction of heavier 
computations. 

4.3 Compatibility conditions 

The knowledge of the angular momentum vector and of the energy at a given 
time allows us to compute the Keplerian elements 

a, e,I,Q . 

In fact the semimajor axis a and the eccentricity e can be computed from the 
energy and the size of the angular momentum through the relations 

£ = ~, ||c||=£V«(l-e 2 ); 

the longitude of the node f2 and the inclination / are obtained from the direction 
of the angular momentum 

c = (sin ft sin /, — cos f2 sin /, cos J) . 

The two attributables Ai,A 2 at epochs ii,t 2 give 8 scalar data, thus the problem 
is over-determined. From a non-spurious pair p 2 ), solution of we obtain 
the same values of a, e, /, fi at both times ti, i — 1,2, but we must check that 
the orbit is indeed the same, that is check the compatibility conditions 

wi=wa, l 1 =t 2 +n{i 1 -i i ), (16) 

where uj\ , u 2 and t\ , i 2 are the arguments of perihelion and the mean anomalies 
of the body at times t\,t 2 and n — fca~ 3 / 2 is the mean motion, which is the 
same for the two orbits. The first of conditions (|16|) corresponds to the use 
of the fifth integral of the Kepler problem, related to Lenz-Laplace's integral 
vector 

T 1 • r 

L= p rxc -R' 

Indeed the compatibility conditions (|TB)) can not be exactly satisfied, due to 
both the errors in the observations and to the planetary perturbations. Actu- 
ally the latter are important only when the observed body undergoes a close 
approach to some planet in the interval between t\ and t 2 . Thus we may be 
able to discard some solutions, for which the compatibility conditions are largely 
violated. Nevertheless, we need a criterion to assess whether smaller discrep- 
ancies from the exact conditions (|16[) are due to the measurement uncertainty 
or rather due to the fact that the two attributables do not belong to the same 
physical object. This will be introduced in the next section. 
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5 Covariance of the solutions 



Given a pair of attributables A = (Ai,A2) with covariance matrices ^An^A 2 j 
we call R = (pi, pi, p 2 , P2) one 01 the solutions of the equation *(R;A) = 0, 
with 

*(R;A)=f yi-VfZ^l^) . (17) 
V £ -t(PuPi) -Zl(P'2,P2) ) 

We can repeat what follows for each solution of <&(R; A) = 0. 

Let R = R(A) = (fti(A),ft 2 (A)), where Hi(A) = (pi(A), Pi(A)) for i = 
1,2. If both the elements (Ai,1Zi(A)), (A2,H2(A)) give negative two-body 
energy orbits, then we can compute the corresponding Keplerian elements at 
times 

i i = i i (A)=U-^-, * = 1,2 

c 

through the transformation 

(a,8,a,6,p,p) = (A,1Z) n> £ K e P {A,Tl) = (a,e,I,Cl,ui,£) . 

We have, for example, a smooth function 

cjj = LJi(A) = u(Ai,TZi(A)) , i = 1,2 

and similar functional relations for a, e, /, 0, £ Actually by construction, we 
have ai — 0,2, ei — e2,h = h,^i — ^2- we denote by a the common value of 
ai and 02. We use the vector differences 

Ai, 2 = (Aw, A*), 

where Au> is the difference of the two angles u>i and 0J2, A£ is the difference 
of the two angles li and £ 2 + n(fi — £2) and n = fca~ 3 / 2 is the mean motion 
of both orbits. Here we compute the difference of two angles in such a way 
that it is a smooth function near a vanishing point; for example we define 
Auj = [wi — L02 +7r(mod 2ir)] —tt. With this caution, the vector Ai ; 2 = Ai j2 (A) 
represents the discrepancy in perihelion argument and mean anomaly of the two 
orbits, comparing the anomalies at the same time ti. We introduce the map 

* : ([-tt.tt) x xl 2 ) 2 ^ [_ 7r , 7r ) x (-J, JxR^l+xlx^xS 1 

(Ai,A 2 ) = A^ ¥(A) = {Ai,Hi,A 1>2 ) , 

giving the orbit (Ai , IZi (A)) in attributable elements at time ti (the epoch of the 
first attributable corrected by aberration), together with the difference A 12 (A) 
in the angular elements, which are not constrained by the angular momentum 
and the energy integrals. By the covariance propagation rule we have 



(18) 
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where 



9A 



I 







dTZi 


dTZi 




dAi 


dA 2 


and 


0A h2 


9Ai, 2 




[ dAi 


8A 2 J 





r Al o 
o r A2 



so that the covariance of '3/ (A) is given by the 8x8 matrix 

r 



*(A) 



Ai 



1 TZi J-Ki.Ai^ 

rAi i2 ,^i rA 1>2 ,Ki rA 12 



where 



FAi,TIi — ^Ai 



dTZi 



1 T 



dAi 



LAi,A 1>2 — ^Ai 



r - dUl v 



dA h2 



dAi 



dA 2 



r^ 2 



_ dAi 

aAi, 2 

dA 2 



1 T 



— J- Ai,Ki J iAi.j,^ -l^Ai ji ^ A 12 ,Ri Ri ; A 1|2 



and 



r - Ml r 



Mi 
<9A 



r - dUl v 



-i T 



OA 



9A li2 



OA 



The matrices f? 1 , «, j = 1, 2, can be computed from the relation 



dA 



d —(A) 



<9$ 

M (R(A),A) 



ax( R ( A )' A ) 



We also have 

— = ^AMA)) + -(^(A))-^ _ _ ( ^ (A)) _Li 
<9Aw <9w /4 /AXX 97ei(A) aw,, <9w /4 /ANN <97e 2 (A) 



14 



and 



9AI 



dAl 
dAl 



^(A^ l( A)) + -(A )Kl (A))^-i 

3 n 
2a 
n 

c 



da , , m , . da , , ^ . A .. 97^i (A) 



[ti(A)-t 2 (A)] 



9p 



Mi 



L (A) 



9p 2 



(A) 



Mi 

^ M 7? fA ^i( A ) 

3 n 
2a 



^ , „ 9^ . „ . „97e 2 (A) 

^(^^ 2 (A))--(^^ 2 (A))^-i 



^(^ 2 (A)) + -(^ 2 (A))^ 



[fi(A)-t 2 (A)] 



9p 



M- 



L (A) 



9p 2 
Mi 



(A) 



5.1 Identification of attributables 

One important step is to decide if trying to link Ai, A2 has produced at least 
one reliable orbit, so that we can state the two sets of observations defining the 
A4, i — 1,2 may belong to one and the same solar system body. Neglecting the 
unavoidable errors in the observations and the approximations made both with 
the interpolation to compute Ai, A2 and with the use of a two-body model, 
if the observations belong to the same solar system body, then Ai i2 (A) = 0. 
We need to check whether the failure of this condition is within the acceptable 
range of values which is statistically expected to be generated by the errors in 
the available observations. 

The marginal covariance matrix of the compatibility conditions is 



r Al . 



9A 1|S 
9A 



9Ai,i 
9A 



The inverse matrix C Al - 2 = r A * defines a norm || • ||* in the (Aw, At) plane, 
allowing to test an identification between the attributables Ai , A2 ■ the test is 



|Ai >2 ||2=Ai, 2 C7^A?; 2 < Xr 2 



(19) 



where Xmax is a control parameter. The value of the control could be selected on 
the basis of x 2 tables, if we could assume that the observations errors are Gaus- 
sian and their standard deviations, mean values and correlations were known. 
Since this hypothesis is not satisfied in practice, the control value Xmax needs to 
be selected on the basis of large scale tests. Note that, for each pair of attributa- 
bles, more than one preliminary orbit computed with the method of Section |4] 
could pass the control (|19p ; thus we can have alternative preliminary orbits. 
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5.2 Uncertainty of the orbits 



The methods explained in Section @] also allow to assign an uncertainty to the 
preliminary orbits that we compute from the two attributables. A solution 
(•At, 72-1 (A)), in attributable elements, has the marginal covariance matrix 



The preliminary orbits obtained by the other available algorithms do not pro- 
duce a nondegenerate covariance matrix: this is usually computed in the differ- 
ential correction step of the orbit determination procedure. With the algorithm 
of [13] a covariance matrix may be defined, but it is not positive definite. The 
advantage of having a covariance matrix already from the preliminary orbit step 
could be important in two ways. First, the covariance matrix describes a con- 
fidence ellipsoid where a two-body orbit, compatible with the observations and 
their errors, can be found. The size of this ellipsoid can provide useful hints 
on the difficulty of the differential corrections procedure. Second, even if the 
differential corrections are divergent, the covariance matrix of the preliminary 
orbit can be used to compute a prediction with confidence region, allowing for 
a planned recovery, for assessment of impact risk, and so on. 

6 A test case 

We show a test of the linkage procedure using the attributables 

At = (0.2872656,0.1106342,-0.00375115,-0.00167695), 
A 2 = (0.2820817,0.1086542,0.00514465,0.00215975) 

of the asteroid (101878) 1999 NR 23 at epochs it = 54000, i 2 = 54109 respectively 
(time in MJD). The values of the components of A\,A2 are in radians/radians 
per day. They have been computed from two groups of observations, separated 
by more than 100 days, made from two different observatories: Mauna Kea 
(568) and Mt. Lemmon Survey (G96). From the known nominal orbit of this 
asteroid we obtain the values 



of the topocentric distance (in AU) at the two mean epochs of the observations. 

In Figure[T]we show the intersections between the curves defined by p(p\ , p 2 ) 
and q(pt,p2). By solving the corresponding problem @ with the method de- 
scribed in Subsection 14.11 we find the 6 positive pairs of solutions (pi,p 2 ) dis- 
played in Table [TJ 

After removing solution 1 (with both components very small), the spurious 
solution 6 (not satisfying (JT))) and the spurious solutions 3 and 5 (not satisfying 
([6])), we are left with the values labeled 2 and 4 in Table [TJ Note that, even if 
solutions 2 and 3 look close, they are far apart enough to select only one of them 




Pi = 1.0419 



= 2.0485 
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(101878) 1999 NR23 , A = 109 d 




0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 

Pi 



Figure 1: Intersections of the curves p = 0, q = (solid and dashed, respectively) 
in the plane pi, pi for the asteroid (101878) 1999 NR23: the asterisk corresponds 
to the true orbit. 

as a good solution. For the remaining solutions 2 and 4 we succeed in computing 
Keplerian orbits, that we list in Table [2] The values of a, e, /, are the same 
for each pair: this is due to the fact that each pair of orbits shares the same 
angular momentum and the same energy. The value of the identification norm 
HAi^H*, also shown in Tabled strongly suggests to select the second solution. 
The results appear pretty good, in fact the differences with the true solution 
are of the order of 3 x 10 -5 AU and the errors in the Keplerian elements are 
comparable with the planetary perturbations; this is intrinsically bound to the 
use of a two-body approximation. 

7 Numerical experiments with simulated obser- 
vations 

We have tested our identification method with the DFT algorithm, explained 
in Subsection 14. 1[ using simulated observations of objects in a solar system 
model. The data have been given to us by R. Jedicke and L. Dennau from the 
Institute of Astronomy, University of Hawaii, and the data quality resemble the 
one which should be achieved by the Pan-STARRS telescope when it will be 
fully operative. The RMS of the observations vary from 0.01 to 0.02 arcsec, 
that is rather optimistic for the current surveys. The current astrometric data 
quality of the Pan-STARRS 1 telescope is such that the RMS of the residuals for 
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Pi 


P2 


1 


0.0059 


0.0097 


2 


0.7130 


1.4100 


3 


0.7045 


1.3933 


4 


1.0409 


2.0517 


5 


1.1659 


2.2952 


6 


1.4246 


2.7968 



Table 1: Solutions of the system © for (101878) 1999 NR 23 . 





1 


2 


l|Ai, 2 ||* 


487.65806 


0.19505 


a 


6.87384 


6.87384 


2.25828 


2.25828 


e 


0.81798 


0.81798 


0.19787 


0.19787 


I 


0.51733 


0.51733 


0.59995 


0.59995 


n 


156.55521 


156.55521 


156.42531 


156.42531 




144.68146 


321.78289 


144.39580 


145.26330 


i 


4.66178 


355.27766 


47.75173 


78.65378 


t (MJD) 


53999.8205 


54109.1368 


53999.8186 


54109.1331 



Table 2: Keplerian elements (angles in degrees) corresponding to the pairs 
[PXiPt) labeled with 2 and 4 in Tabled] The value of || Ai^H* is shown for 
the two solutions. 

well determined asteroid orbits is between 0.11 and 0.13 arcsec. Better results 
should be achieved when the astrometric reduction of asteroid detections will 
be performed with respect to a catalogue generated by the Pan-STARRS survey 
itself. 

These observations cover 31 observing nights, in three consecutive lunations and 
are grouped into tracklets. Each tracklet is composed by observations presum- 
ably belonging to the same object and covering a short arc: some of them are 
false (e.g. join observations of different objects). From each tracklet we can 
compute an attributable. We have first applied to the database of tracklets the 
identification procedures defined in [T3], [7]. Then we have tested our method 
on the leftover database, for which the previous procedure has failed. These 
remaining observations corresponds to 19441 objects, and 24590 tracklets, but 
only 4132 objects have at least two tracklets, that is a necessary requirement for 
the application of our method. The hyperbolic orbits have been removed from 
the solar system model: in fact our current method does not search for them, 
but it could be easily modified to include their orbit determination. To each 
accepted preliminary orbit obtained from a pair of attributables we apply the 
differential corrections, using all the observations at our disposal, to compute a 



18 



least squares orbit with its covariance matrix^ 

To reduce the computational complexity, we need to define a filter selecting 
the pairs of attributables which we try to link. In Section 17.11 we describe the 
two filters we have used in processing the simulated data. 

7.1 Filtering pairs of attributables 

7.1.1 First filter: guessing the second angular position. 

A first simple way to discard pairs of attributables at epochs t\ , ti is to constrain 
the time span St = ?2 — it- we require 

Stmin < St < 5t max (20) 

for suitable positive constants St m i n ,St max - In our experiment we have used 
Stmax = 99 days and St m i n = 0.5 days, that practically means we have tried to 
link attributables obtained in different nights. For each given pair of attributa- 
bles at epochs t\ , i% fulfilling (|20|) we consider for i = 1, 2 the corresponding 
proper motions m and the mobile bases {pi,Vi,&i}, defined in Section [21 We 
want to use one of the proper motions, say r\\, to bound the region in the sky 
where we could recover the object at the other time i% . 

Let us form the orthogonal matrices V\ = [pi|vi|ni] and Vi — [p2|v 2 |n 2 ]: 
these are rotation matrices to the mobile bases {pi,Vi,hi}, i = 1,2. Let 
denote the rotation of an angle around the unit vector e. Then R m st fh = 
Vi R m sti V\ (z is the third unit vector of the reference frame defining our 
rectangular coordinates) is the parallel transport matrix along the geodesic on 
the unit sphere defined by Ai to time i^] hence P12 = -fi^stniPi is the predicted 
observation direction at time t 2 , assuming the trajectory is a great circle and 
the proper motion is constant. By exchanging the order of the two attributables 
we can compute R- V2 sta 2 — V2 R-^stz V 2 and P21 = R-r] 2 5th 2 P2, that is the 
prediction at time t\. We use the metric 

d{pi,fa) = min{pi2,p 2 ,P2i,Pi} , 

that is the minimum between the two angular differences, discarding pairs of 
attributables that give rise to a large value of this metric. 

Note that the proper motion does not vary too much in the time interval 
between the two attributables provided St max is small enough; thus, if we want 
to use large values of St, we have also to allow large values of the metric d. 

7.1.2 Second filter: symmetric LLS fit. 

Given the two attributables Ai , Ai at times t\ , ti we perform a quadratic ap- 
proximation of the apparent motion on the celestial sphere S 2 by using a Lin- 
ear Least Squares (LLS) fit. The apparent motion is given by the functions 

5 We use the preliminary orbit at time t\ as starting guess for the differential corrections. 
We could also use the orbit at time ti, or an 'average orbit' at time (ti + ?2)/2. 
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a(t), 6(t). We approximate a(t), S(t) with second degree polynomials whose co- 
efficients are derived from a least squares fit. We denote the approximating 
quadratic functions as 



a(t) = a q + a q (t-t) + ^a q (t-t) 2 . 



5(t)=5 q + 5 q (t-t) + ^6 q (t-t) 2 , (21) 



where t = |(ix + t 2 ) is the mean of the times of the attributables. The corre- 
sponding time derivatives are 

a(t) = a q + a q (t - t) , S(t) = 6 q + S q (t - t) . 

We want to determine the 6 quantities a q , a q , a q , 5 q , S q , 5 q using the data coming 
from the attributables. The vector of residuals is 



£ = (A 1 -A(t 1 ),A 2 -A(i 2 )y , 
with A(t) = (a(t),5(t),d(t),S(t)). 

Given the covariance matrices r^ 1 , T_^ 2 associated to the attributables we use 
them to weight the residuals in the definition of the target function: 



Q($)=y-w$., 

We introduce the notation 



where W = 



T M 



x = (a g ,ag,a g ,S q ,S q ,5 q ) T , A = (A ll A 2 ) T ■ 

The value of £ = £(x) that minimizes the target function is obtained by solving 
the normal equation 



Cx = —B WX , 



and the matrix B has the form 



where B = 



<9x' 



C = B T WB . 







B 



B 2 



with Bi 



( 1 (U-t) \{U-t) 2 

1 {U-t) \{U-t) 

1 {U-t) 

\ 1 {U-t) J 



for i = 1,2. Once the value of x is given, we compute the residuals £(x) and 
use the norm y/Q(£) to decide which are the pairs of attributables (^1,^2) to 
discard. We also discard the pairs giving rise to a large value of the quantity 



K qV q 



1 r 

'hi 



{5 q a q - a q 5 q ) cos5 q + a q (rj q + S q ) smd q 



where rj q = y 5^ + a 2 q cos 2 S q and n q is the geodesic curvature (see [TB] , Chapter 
9). 
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7.2 Results 



The accuracy of the linkage method can be measured by the number of true 
identifications over the total number of identifications found. This computation 
includes duplications due to alternative solutions. The total number is 3906 
and the true ones (that may be related to the same object if it has more than 
2 tracklets) are 3144, i.e. 80.5% of the total. We could eliminate almost half 
of the 762 false identifications by lowering from 0.15 to 0.06 arcsec the control 
on the RMS for acceptable orbits after differential corrections: but this would 
make us lose 102 true identifications. 

In Table [3] we show the efficiency of the linkage procedure, that is we write 
the number of objects for which at least a pair of tracklets has been correctly 
linked, giving the details for the MB (Main Belt) and the NEO (Near Earth 
Object) class. As expected, the efficiency appears greater if there are three 
tracklets that can be pairwise linked. We stress that we have tested our method 



with 2 tracklets in 2 nights 


Total 


Found 


Lost 


all 
MB 
NEO 


1074 
1038 
19 


963 
947 
9 


89.7% 
91.2% 
47.4% 


111 
91 | 
10 


10.3% 
8.8% 
52.6% 










with 3 tracklets in 3 nights 


Total 


Found 


Lost 


All 

MB 
NEO 


214 
197 

3 


205 
196 
2 


95.8% 
99.5% 
66.7% 


9 
1 
1 


4.2% 
0.5% 
33.3% 



Table 3: Efficiency of the identification procedure. 

with data for which the other available methods in [13], [7] could not perform 
the linkage. 

8 Conclusions and future work 

We have investigated an orbit determination method that is based on the inte- 
grals of the Kepler problem and is suitable to be used with modern data sets 
of observations. With the recent technologies we can collect a large number of 
tracklets in each observing night and it is difficult even to relate two tracklets 
of different nights as belonging to the same observed object. We have defined 
a linkage procedure between tracklets, allowing to compute preliminary orbits 
with covariance matrices. An interesting feature that comes out from our nu- 
merical experiments is that this method appears to work also when the time 
span between the two attributables is large, hence it can be used in cases where 
other linkage methods fail. The efficiency and performance of the algorithm 
explained in Subsection 14.11 have been studied with a large scale test. Therefore 
this method can be important for two kinds of applications: 1) to recover ob- 
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jects whose orbit could not be computed with either the classical or the modern 
known algorithms; 2) to design the scheduler of new surveys planning a smaller 
number of observations for each object. 

The number of alternative solutions of the problem deserves a deeper inves- 
tigation, however we expect that the acceptable ones should often be much less 
than 48, the total degree of the polynomial system ©. Moreover the perfor- 
mance of the second algorithm to solve ([9]), described in Subsection l4.2[ has not 
been tested yet: we would like to perform further experiments to decide if it 
allows to decrease the computation time. 
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